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Abstract 

We investigate the spatial and temporal features of dense contaminant plumes dynamics in porous materials. Our analysis is 
supported by novel experimental results concerning pollutant concentration profiles inside a vertical column setup. We describe 
the experimental methods and elucidate the salient outcomes of the measurements, with focus on miscible fluids in homogeneous 
saturated media. By resorting to a finite elements approach, we numerically solve the equations that rule the pollutants migration 
and compare the simulation results with the experimental data. Finally, we qualitatively explore the interfacial dynamics behavior 
between the dense contaminant plume and the lighter resident fluid that saturates the column. 

Key words: Variable-density contaminant dynamics. Homogeneous porous media, Non-Fickian transport 
47.56.+r, 47.20.Ma 
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1. Introduction 

The migration of dense pollutant plumes through porous 
materials is key to mastering such environmental and tech- 
nological challenges as site remediation, safety assessment 
tor waste reposi tor ies and enhanced pet ro leum recov - 
ery jOphori 2004 lYang and Edwards 2000t iHolm 19861) . 
'to name a few. Yet, despite be i ng a well-known and lon g 
studied problem (^ Wooding 1962 : Bachmat and Elrick X970'), 



>^ it kee ps raising r nany conceptual a s well as pra ctical is 



sues jBear 19721: lOe Marsilv 19861: ISahimi 19951): for an 



overview of recent advan ces, see, e.g., fSim mons et al. 2001 



iDiersch and Kolditz 20021) . When contaminants are suffi- 
ciently diluted, flow and concentration fields do not affect 
each other, so that pollutants dynamics can be addressed by 
first determining the pressure distribution in the traversed 
region (which imposes velocity within the pore network) 
and then separately computing the concentration profile 
due to advection and dispersion mechanisms. Under such 
conditions, transport is Fickian, i.e., r uled by th e standard 
advection-dispersion equation (ADE) (iBear 1972h . at least 
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for homogeneous media. On the contrary, when contaminant 
plumes are sufficiently concentrated the pressure field will 
depend on the fuid density, which in turn is afl'ected by the 
contaminant concentration. Nonlinear effects come thus into 
play, by virtue of th e couplin g between flow and transport (see , 
e.g.. jSchincariol an d Schw artz 1990l: IWeltv and Gelhar 1991 



Weltv and Gelhar 1992i; iWooding 1969h "). Based on the den 



sity (and in general also viscosity djiao and Hotzl 2004)) 
contrasts between the pollutants-loaded fluid and the 
resident fluid that permeates the porous material, insta- 
bilities (fingerings) may appear at their mutual inter- 



face dManickam and Homsv 19951: IWooding et al. 1997 



Wooding 19691) . The relevance of such instabilities to the 
overall contaminant dynamics is determin ed by the inten 
sity of density an d viscosi ty contrasts (iBacri et al. 1991 



Jiao and Hotzl 2004^ 

Welty and Gelhar 1991j 

Oltean and Bues 2002l:]Rogerson and Meiburg 1993h . 



Schincariol a nd Schwartz 199C; 

2008; 



jDalziel et al. 



High densities and/or density gradients are encountered 
when either the contaminant itself is strongly concentrated 
at the source, or the plur ne flows through regions that 
are rich in salt (brines) ( Hassaniz adeh and Leijnse 19951; 
Schincariol and Schwartz 199o[ ISchotting et al. \99% . 
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Figure 1: Scheme of the BEETI experimental setup for downwards and up- 
wards tracer injection. 



In particular, the latter case might become a ma- 
jor con cern for radioactive waste disposal near salt 
domes ( Hassanizadeh and Leiinse 1995h . Variable-density 



migration affects transport by altering both the time scale 
(average disp lacement) and the spatial extent (spread) of pollu - 



tants plumes (ISimmons et al. 2001l : lDiersch and Kolditz 20021) . 



Therefore, in such conditions the plume dynamics can not 
be described by resorting to the stand ard APE model, which 



does not account for density effects (iWood et al. 20041) . An 
exact quantification of the pollutants average displacement 
and spread depends on correctly capturing the underlying 
physics of the above mentioned processes: several studies have 
been performed to this aim, covering homogeneous saturated 
and heterogenous unsaturated materials, density as well as 
viscosity contrasts, and dif ferent geometries and column 



Wood etal. 2004; 


Jiao and Hotzl 20041 Dalziel et al. 2008: 


Oltean and Bues 2002| 


D'Anaeloet al. 2008; 


Simmons et al. 20021 


TchleDi et al. 1993" iDe Wit et al. 2005; 


Zimmerman and Homsv 1992; Ruith and Meiburs 2000 ; 


Camhi et al. 2000t 


Debaca et al. 200 ll Seon et al. 2004 ; 


Debacq etal. 2003) 




A central outcome is that density 



effects can play a major role in affecting the fate of the 
contaminant plum e even in homoge neo us and saturat ed 



porous media (iJiao and Hotzl 2004t Wood et al. 2004 



Simmons et al. 2002r lOltean and Bues 2002h . Moreover, 
it has been shown that even modest density differences 
with respect to the reference fluid (which is frequently 
freshwater) may lead to m easurable deviations from the 
usual Fickian behavior ( Schincariol and Schwartz 19901; 



i. 



Hassanizadeh and Leiinse 19951) . 

The aim of this work is to explore the spatial and tempo- 
ral features of dense pollutants transport in porous media, with 
focus on the case of miscible fluids, i.e., fluids having a sin- 
gle phase but different solute concentrations. Our investiga- 
tion is supported by novel experimental measurements of con- 



taminant profiles inside homogeneous saturated columns, ob- 
tained by means of X-ray spectroscopy. As detailed in the 
following. X-rays measures allow for two complementary in- 
formations: on one hand, the contaminant profiles as a func- 
tion of time, at a given location; on the other, the spatial evo- 
lution of such profiles along the column. We are thus able 
to fully characterize contaminant dynamics by experimental 
means, which is especially important when the concentration 
profiles are time-dependent. The effects of unstable interfa- 
cial dynamics on the spatial and temporal evolution of con- 
taminant profiles have received so far only limited attention, 
whereas studies are usually focused on the mixing proper- 
ties at the inter face between two laye rs of semi- infinite exten- 
sion (see, e.g., dWood et al. 20 04; D e Wit et al. 2 005) and ref- 
erences therein). In groundwater contamination, the spill ex- 
tent is often finite, since the source is limited in space and 
time (De Wit et al. 2005) ; based on this observation, our anal- 
ysis will focus on finite-duration contaminant injection. In this 
case, dynamics is transient: stable and unstable fronts (between 
injected and resident fluid) will in general be present at the same 
time. 

This paper is organized as follows: in Sec.|2]we describe the 
experimental setup and outline the salient features of the con- 
taminant concentration profiles along the column. Then, we 
review the underlying physical equations in Sec. [3] and address 
their numerical integration by resorting to finite elements mod- 
eling in Sec.S] Simulation results are successively compared to 
experimental measurements in Sec. |5] The qualitative behav- 
ior of the interfacial instabilities and the interplay of the differ- 
ent components governing the physical system are explored in 
Sec.|6] Conclusions are finally drawn in Sec.]?] 

2. Experimental methods and results 

2.1. The experimental device 

The BEETI experimental setup has been specifically con- 
ceived to assess the spatial and temporal dynamics of tracers 
in relation to the geometry and the physical-chemical condi- 
tions of porous media. The system is composed of a vertical 
polycarbonate column, a controlled hydraulic circuit that en- 
ables tracers migration within the column, and a scanning X- 
ray spectrometer for measurements (Fig. [Hi- The column has 
height // = 80 cm and internal radius r = 2.5 cm, so that the 
aspect ratio is Hjlr - 16 » 1. This experimental setup allows 
for downwards as well as upwards fluid injection, and several 
kinds of flow regimes and porous materials can be tested, at 
various saturation and/or heterogeneity conditions. 

The scanning X-ray spectrometry system (SXSS) is based 
on an X-ray generator with a tungsten source. The emerging 
beam is filtrated by a neodymium window that selects two en- 
ergy ranges (20 - 40 keV and 50 - 75 keV, respectively). After 
traversing the sand column and its water and tracers content, 
the beam is measured by a Nal detector Both source and detec- 
tor are kept in position by a controlled rack rail that displaces 
the SXSS along the column. At specified spatial locations, the 
SXSS waits 60 seconds (the counting time) and measures the 
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Figure 2: Contaminant concentration curves Cc(f) at tlie last measure point 
along the column, C = 11 cm, as a function of time [h] . Downwards injection at 
increasing molarities C""'' = 0.05 (solid line), 0.1 (crosses), 0.2 (triangles) and 
0.5 (circles) mol/L. 



transmitted dichromatic X-rays. The X-ray countings received 
by the detector depend on the thickness of the traversed layer, 
the nature of traversed phase (fluid, porous material or tracer) 
and the counting time. It is possible to discriminate the diflferent 
components within each column layer by resorting to the Beer- 
Lambert's law, which is used to convert the transmitted beam 
intensity to physical quantities (such as bulk density, porosity, 
or tracer concentration): 



R = Ro exp 



■Z 



(1) 



Here Rq and R denote the number of photons emitted and de- 
tected, respectively, per unit time; / is the index of each phase 
contained in the medium (sand, water, solute, ...); yu, [m^/Kg] 
is the mass attenuation coefficient, which depends on the phase 
composition and on the X-ray energy; p, [Kg/m^^*] is the den- 
sity, and Xi [m] the thickness of the phase. The parameters yu,- 
and Rq are determined by calibrating the SXSS with specific 
materials specimens, with perfectly measurable thickness and 
X-ray attenuation very close to that of water or sand. One of the 
specimens is used as reference for the calibration of Rq, so to 
avoid saturation effects of the Nal detector, due to strong source 
intensity. Tracer concentration is determined by rewriting the 
Beer-Lambert's law as a function of the molecular attenuation 
coefficient, which is linearly proportional to the tracer concen- 
tration solution. The X-ray transmitted countings allow quan- 
titatively accessing the tracer concentration inside the column 
(as a function of time), at a distance { from the inlet: we denote 
this quantity by C({t). Given the spatial resolution of the X-ray 
device, Ce(t) actually corresponds to a spatially-averaged mea- 
sure over the transverse direction (i.e., the column section), for 
any given height {. 

In order to obtain a homogeneously packed porous medium, 
the dry sand was poured continuously and packed while mak- 
ing the column vibrate. Solution saturation was imposed start- 



Figure 3: Contaminant concentration curves Cf(f) at the last measure point 
along the column ( = 11 cm, as a function of time [h]. Upwards injection at 
increasing molarities C"'°' = 0.05 (solid line), 0.1 (crosses), 0.2 (triangles) and 
0.5 (circles) mol/L. 



Molarity [mol/L] 


Density [Kg/m^^] 


Ap/po [-] 


0.05 


1004.3 


0.006 


0.1 


1010.4 


0.012 


0.2 


1022.4 


0.024 


0.5 


1058.5 


0.060 



Table 1 : Density vaiiations for KI in water, at increasing solution molality. 



ing from the column bottom. Medium homogeneity and satu- 
ration were assessed by dichromatic X-ray spectrometry, which 
allows estimating bulk density and porosity along the column 
with a resolution of 5 mm in the longitudinal direction. 

Upwards experiments (Fig. [1] left) are performed by inject- 
ing background and tracer solutions from the lower entrance. 
Conversely, the fluids are injected from the upper entrance 
for downwards experiments (Fig. [T] right). The solutions are 
drained by a pump which regulates the solution flow rate and 
the tracer injection at the column inlet. The steady state Darcy 
flow q is also verified by weighing the outgoing solution. Elec- 
tric conductivity cells are adopted on-line at the column in- 
let and outlet, so to monitor tracer displacement through the 
column. At the outlet of the column, C(=H{t) coincides with 
the breakthrough curve, which is the most frequently measured 
quantity in contaminant migration experiments and provides in- 
formation about the tracer mean displacement over the entire 
column length. 

2.2. Experimental conditions 

In the following, we refer to a saturated column filled with 
homogeneously mixed Fontainebleau quartz sand, with bulk 
density 1.77 + 0.01 g/cnr' and average grain diameter 200 fim. 
The sand was previously washed in pure water (MilliQ) and 
equilibrated with a KCl 10"^ M solution, so to minimize the 
ionic exchange between solid sites surface and tracers. 
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Figure 4: Downwards injection at a reference molarity C""' = 0.2 
mol/L. Contaminant concentration curves Cf(f) measured at sections ( = 
7.7,23.1,38.5,46.2, and 77 cm (from left to right), as a function of time [h]. 
Lines have been added to guide the eye. 



Figure 5: Upwards injection at a reference molarity C"'°' = 0.2 mol/L. Contam- 
inant concentration curves Ce(t) measured at sections f = 7.7,23.1,38.5,46.2, 
and 77 cm (from left to right), as a function of time [h]. Lines have been added 
to guide the eye. 



The average porosity is = 0.333 + 0.05, as measured by the 
SXSS (to be compared with the value = 0.328 as obtained by 
weighing the sand poured into the column). The dispersivity is 
a = 0. 1 cm, as from previous experimental runs (at higher flow 
rates) on the BEETI device. All measurements are performed at 
constant room temperature T = 20°C. The reference saturating 
solution contains KCl at a molar concentration of C™"' - 10^^ 
mol/L (molar mass equal to 74.5 g/mol), so that the correspond- 
ing reference density is po - 998.3 Kg/m^ at T = 20''C. The in- 
jected tracer is KI (molar mass equal to 166 g/mol), at various 
molar concentrations. The examined experimental conditions 
are presented in Tab.[Tl together with the corresponding injected 
fluid densities p anddifferential densities Ap/po - (p - Pa) I Pa, 
as derived from (lOtaka 1985D . 

The stationary flow \s q = 2 cm/h and the duration of the 
flux step injection is T = 3 h. In most experimental runs, the 
measure points are located at a distance C - 1.1, 23.1, 38.5, 
46.2 and 77 cm from the column inlet. 

The experimental conditions are such that clogging or forma- 
tion of colloidal particles, which could alter the interpretation 
of the obtained results, can be excluded. Chemical reactions or 
sorption/desorption phenomena can be ruled out as well. 

2.3. Experimental results 

In the following, we provide some representative results 
that can best illustrate the set of measures performed with the 
BEETI device. In Figs. |2] and [3] we start by displaying the 
concentration profiles corresponding to the last measure point 
along the column as a function of time, for downwards and up- 
wards injection, respectively. Here and throughout the text, 
concentration profiles have been normalized to their respec- 
tive areas. Curves are shown for increasing values of the con- 
centration molarity, ranging from C"'"' = 0.05 to C'"°' - 0.5 
mol/L. At weak molarity, the standard Fickian (symmetrical) 



profile is recovered. As C™"' increases, the curves become in- 
creasingly skewed and appreciably deviate from the Gaussian 
shape. The sign of the skewness depends on the flow direc- 
tion: downwards injection leads to negatively skewed profiles 
(i.e., contaminant flows out earlier than expected for a Fick- 
ian flux: see Fig. |2]i, whereas upwards injection leads to pos- 
itively skewed profiles (i.e., contaminant flows out later than 
expected: see Fig.[3]l. These eff'ects on transport can be mainly 
attributed to gravity (via the density coupling). An unstable 
front appears where a dense fluid overlies a lighter one, and 
the macroscopic effect is an enh anced difFusivity (i.e., spread) 
at the interface (I Wooding 19691): this mechanism is known as 
Rayleigh-Taylor instability (Tavl or 1950l) . Actually, the inter- 
play between mechanical dispersion due to the structure of the 
porous medium and the additional dispersion due to interfa- 
cial i nstabilities is ke y to understanding variable-density trans- 
port dWooding 1969b . 

Then, in Figs.|4]and|5]we display the evolution of concentra- 
tion profiles measured at various heights along the experimental 
device, for downwards and upwards injection, respectively. In 
fact. X-rays spectrometry allows accessing contaminant distri- 
bution inside the column. Profiles are shown at a fixed mo- 
larity C"'"' = 0.2 mol/L, which provides a representative case. 
From these results it appears that the skewness of the curves 
increases along the column, in the direction dictated by the un- 
stable front. In other words, the portion of the plume that moves 
faster or slower than the bulk progressively increases. This is 
because (for the experimental conditions considered here) in- 
stabilities can grow, overcoming the smoothing effects of the 
dispersion induced by the porous material. We remark also that 
the concentration profiles are almost symmetric upon reversing 
flow direction, which underlines the prominent role of gravity 
in determining the pollutants dynamics. 

Finally, an important question concerns the reproducibility of 
the experimental results presented here. As a general remark. 
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Figure 6: Contaminant concentration curves Cc(f) at tlie last measure point 
along the column, C = 11 cm, as a function of time [h]. Downwards injection 
at C'""' = 0.5 mol/L, for 5 experimental runs. 



Figure 7: Contaminant concentration curves C((t) at the last measure point 
along the column, t = 11 cm, as a function of time [h] . Downwards injection 
at C"'°' = 0.2 mol/L, for 4 experimental runs. 



the presence of instabilities at the interface between resident 
and displacing fluids leads to fluctuations and thus to random- 
ness in the shape of contaminant concentration profiles; this 
holds true even for homogeneous and saturated p orous media, 
such as our column setup (Sim mons et al. 200 lb . Triggering 
of unstable fronts and their interplay with stable fronts will be 
discussed later. At strong molarity (say C'""' = 0.5 mol/L or 
stronger), the effects of the unstable front dominate, so that a 
limited reproducibility of the experimental curves is actually 
observed, because of the aforementioned stochasticity. In Fig.|6] 
we display results for the repetition of 5 experimental runs for 
downwards injection at C'"°' = 0.5 mol/L: only a small por- 
tion of the curves is reproducible, and fluctuations are evident. 
Note however that the average displacement and the spread of 
the contaminant plume seem largely conserved, although the 
specific shapes of the profiles vary from run to run. 

At weak molarity (say C™"' = 0.05 mol/L or weaker), con- 
taminant profiles are perfectly reproducible within the limits 
of experimental measures. This is an expected outcome, since 
in this case the flow and transport mechanisms are almost de- 
coupled and the plume dynamics approaches the standard Fick- 
ian behavior Finally, at intermediate molarity, we experimen- 
tally observe a gradual transition between these two regimes. 
In Fig. |7]we display the repetition of 4 experimental runs for 
downwards injection at C'""' = 0.2 mol/L: a good reproducibil- 
ity of the concentration profiles is found, as the curves are 
almost superposed, although fluctuations are already visible. 
Note in particular that superposition is more remarkable on the 
right portion of the curves, whose shape is determined by the 
stable front. Comparable results have been obtained also for 
the case of upwards injection (not shown here). 

3. The coupled flow-transport model 

We address now the description of the physical model un- 
derlying the observed contaminant dynamics, on the basis of 



the experimental results exposed above. For sake of simplic- 
ity, we commence by representing the actual geometry of the 
vertical column as a flat rectangular {2d) domain Q, of height 
H [cm] and base B [cm]: Q. = [0,B] x [0,//]. We assume 
that this region is homogeneous and isotropic, as being filled 
with well-mixed saturated sand. The dynamics of a contami- 
nant plume injected in Q is ruled by three equations that con- 
dense the physics of the problem, namely fluid mass conser- 
vation, advection-dispersion of the contaminant concentration 
field c(x, t) and Darcy's law for the advection field u(x, t). Con- 
cerning mass conservation, it is customary to introduce the so- 
called Boussinesq approximation, i.e., to assume that the effects 
of density variations in rigid porous media are retained as be- 
ing significant only when appear ing multiplied by the gravity 
constant g (see, e.g., dBear 1972h l. This is indeed the case for 
most examples of dense contaminant transport, where the den- 
sity variations are of the order of a few percents with respect to 



the resident fluid jSimmons et al. 2002 : Simmons et al. 2001 



Dentz et al. 20061) . Under this hypothesis, fluid mass balance 
simplifies to 



V ■ u(x, f) ^ 0, 



(2) 



in the absence of sources and sinks. Mass transport is 
given by the standard advection-dispersion equation, pro- 
vided that the traversed materi al is sufficiently homoge- 
neous (ICortis and Berkowitz 20041) : 



— 0c(x, f) = V ■ [6»D(x, f)V - u(x, f)] c(x, f), 
at 



(3) 



where 6 is the constant average porosity of the medium, u(x, f) 
is the advection field, and the apparent dispersion D includes 
both mechanical dispersion and molecular diffusion. A widely 
adopted expression for D reads (.Bear 197Z) 



D(x, t) 



a|u(x, f)l 



(4) 
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Figure 8: Downwards injection at a reference molarity C 



'Hiol 



0.2 mol/L and 



q = 2 cm/h. Contaminant concentration curves C((f) measured at sections 
t = 7.7, 23.1, 38.5, 46.2, and 77 cm (from left to right), as a function of time 
[h]. Squares represent experimental data, dashed lines model estimates. 



Figure 9: Upwards injection at a reference molarity C"'°' = 0.2 mol/L and 
q = 2 cm/h. Contaminant concentration curves Cc(t) measured at sections 
t = 7.7,23.1,38.5,46.2, and 77 cm (from left to right), as a function of time 
[h]. Squares represent experimental data, dashed lines model estimates. 



where mechanical dispersion is assumed to be proportional to 
the absolute value of the pore velocity u/0 through a constant 
length scale a (the dispersivity) and apparent molecular diffu- 
sion has a constant value Do, incorporating the effects of tor- 
tuosity. In our experimental conditions, it turns out that Do is 
actually negligible with respect to the contribution of mechan- 
ical dispersion: D(x, f) ^ a|u(x, Ol/fi*- Finally , the advec tion 
field is given by the modified Darcy's equation (IBear 1972h 



u(x, t) 



Hc(x, t) 



(5) 



where Vp is the pressure gradient, g is the (vector) gravity con- 
stant (with modulus g), k is the intrinsic permeability of the 
porous material (here assumed to be constant, since the medium 
is homogeneous), //c(x, f) is the fluid viscosity and Pc(x, f) is the 
fluid density. In principle, both density and viscosity depend on 
concentration. Therefore, in order to close the system of equa- 
tions given by |2] |3] and |5] we have to assign the equations of 
state for /ic(x, t) and Pe(x, t). Here, linear constitutive relation- 
ships are adopted, namely 



Pc(x,f) =po[l +ec(x,f)] 



and 



//c(X,f) = /UO [1 +rc(X,0] : 



(6) 



(7) 



which in most cases have been shown to accurately re- 
produce experim entally measured variat io ns (see, e.g., t he 



discussions in (| We Itv and Gelhar 19911: iDentz et al. 2006 



are 



Diersch and Kolditz 2002h ). Ot her functional forms 
also possible (see, e .g.. dHa s sanizadeh and L e ijnse 1995 



Welty and Gelhar 199ll: lOiersch and Kolditz 2002i) and refer- 
ences therein). The values fj.o and po represent the reference 
viscosity and density of the fluid, respectively (when no con- 
taminant is present). The parameters e and y depend on the 



specific injected contaminant species and are defined as e = 
(PQ^)dpc/dc and y = Qify^)dfj.c/dc, evaluated at given exper- 
imental conditions. The magnitude of these parameters ulti- 
mately determines the strength of the nonlinear coupling be- 
tween advection and concentration in Eq.[3] 

Darcy's law |5] can be conveniently rewritten in terms of the 
piezometric head h = p/gpo + y, y being the vertical axis coor- 
dinate: 



u(x, f) = - 



K 



1 + yc{x, t) 



[v/z + ec(x, f)ey] , 



(8) 



where K - kpog/iupi s the hydraulic conductivity of the porous 
medium (IBear 1972h and Cy is the unit vector in the vertical 
direction. Then, combining Eqs. [8] and |2] we finally have 



1 



1 -t- yc(x, t) 



[Vh + ec(x, t)ty] I = 



(9) 



i.e., a single equation for the unknown piezometric head 
h. In the resulting physical system provided by Eqs. |9l [3] 
and |8] the unknowns to be determined are the piezometric 
head h, the advection field u(x, f) and the concentration field 
c{x,t). Some authors mention that the variations of viscos- 
ity are less relevant than those of den sity, for typical transport 
problems in porous media (see, e.g., c We Itv and Gelhar 1991 



Weltv and Gelhar 1 992)): we choose then to set fic - jUo in the 
following. This hypothesis is to be tested a posteriori by com- 
paring model outcomes with experimental measures. 

The system above is complemented by assigning the proper 
boundary and initial conditions, in agreement with the exper- 
imental configuration. On the lateral sides of the column we 
have u • n = and DVc ■ n = 0, n being the normal outwards 
vector This means that there is no flux of matter on the lateral 
sides. Note that this condition also implies V/i • n = on the 
lateral sides. Moreover, we impose a fixed value of the piezo- 
metric head h: it is expedient to set /i = at the outlet of the 
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Figure 10: Downwards injection at a reference molarity C"'°' = 0.2 mol/L 
(cf. Fig.[8). Cumulants Kk[t{i)}, k = 1, 2, of passage times t((') at various column 
heights {. Crosses represent the mean of the passage times (k = 1), circles the 
variance (k = 2) computed from experimental data. Solid (k = I) and dashed 
lines {k = 2) are the model estimates. 



Figure 11: Upwards injection at a reference molarity C™ 



0.2 mol/L 



(cf. Fig.|9). Cumulants a:j [?(/")]. ^ = 1 , 2, of passage times f(f) at various column 
heights £. Crosses represent the mean of the passage times {k = 1), circles the 
variance (k = 2) computed from experimental data. Solid (k = i) and dashed 
lines (k = 2) are the model estimates. 



column. As an initial condition, we assume that the concentra- 
tion field is zero everywhere, c(x, 0) - 0, as the porous material 
is saturated with the reference fluid. At the column inlet, we as- 
sign the ingoing flux in the form of a Robin (mixed) boundary 
condition: 



- [u(x, f) - ^^(x, f)V] c(x, f) ■ n|an„, = ^o(f), 



(10) 



where we have denoted by c?Q,„ the entrance boundary. For the 
experimental runs that we have performed, the imposed flux 
Jo(t) is modulated in time as a finite-duration step injection, 
from t - to r. Finally, we specify a Neumann boundary 
condition, DVc • n = 0, at the outlet of the column. 

4. Numerical solution 

Several alternative approaches have been proposed 
for solving the system of equations |9] |3] and |8] In 
some special cases, it is pos sible to find approxi- 
analytical solutions (lOltean and Bues 2002h . espe- 
for stationary regimes (IDentz et al. 20061) . Aside, 
efforts have been also devoted to the develop- 



mate 
cially 
many 



ment of efi"ective Id mod els (ITardv and Pearson 2006c 



Hassanizadeh and Leiinse 19951: ^[Landman et al. 2007, 

Landman et al. 2007; Egoro v et al. 2005t |Liu and Dane 1996). 



These models are helpful when dealing with column de- 
vices, which can be regarded as being almost- Ic/. In 
general, however, one must resort to numerical solu- 
tions; for instance, Eulerian methods on a discretized 
domain, such as finite differences, finite elements or spec- 
tral methods, either applied to the model above, or to a 
set of an cillary equations for t he fl ow lines (stream func- 



tions) dW eltv and Gelhar 199It IWeltv and Gelhar 1992'; 
Rogerson a nd Meiburg 1993; |Ruith and Meiburg 2000; 

Camhi et al. 20001) . Fully Lagrangian methods are based 



instead on following fluid parcels along their trajectories 
through the porous medium (Zoia et al. 2009). As tra- 
jectories are correlated via the flow-transport coupling, 
this approach becomes computationally expensive for re- 
alistic multi-dimensional domains. Another Lagrangian 
approach. Smoothed Particle s Hydrodynamics (SPH), has 
been recently proposed (iTartakovsky and Meakin 2005 



Tartakovsky et al. 2008), which has shown good potential in 
dealing with Rayleigh-Taylor instabilities in porous media, 
especiall y when stocha sticity is added to describe disper- 
sion (Tart akovsky et al. 2 008). An hybrid approach, combining 
Eulerian and Lagrangian methods, consists in solving the flow 
equations by means of, e.g., finite differences on a grid, and 
then addressing transport via Monte Carlo simulation of con- 
taminant particle s stochastically mov i ng on the grid, at assigne d 



advection field dTchlepi et al. 19931: lAraktingi and Orr 19881) . 



At each time step, the concentration profile derived by particle 
tracking is given as input to the numerical solver, in order to 
update the pressure distribution and re-compute advection. 
Finally, an increasing popular computational tool to tackle 
unstable i nterfacial dyna mics is the Lattice Boltzmann Method: 
see, e.g., dHe et al. 199% and references therein. 

For our aims, we resort to a Eulerian approach and adopt 
the finite elements general-purpose code CAST3M (which has 
been developed at CEA 0) for solving the system |9l [3] and [8] 
above. 

First, the time discretization of the equations is performed via 
a finite differences approach, by introducing a (small) parame- 
ter At as the time step for integration. For sake of simplicity, 
we will consider the backward Euler implicit formula (BDFl); 
the extension to second-order accurate backward difference for- 
mula (BDF2) would be straightforward. In order to cope with 



'CAST3M website: |http://www-cas t3m.cea.fr/cas t3m/index.jsp| 
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Figure 12: Central moment K3[t(£)] of passage times t(£) at various column 
heights {, for downwards (squares) and upwards (triangles) injection at a refer- 
ence molarity C'""' = 0.2 mol/L. Model estimates for the two cases are drawn 
as solid lines. 



the nonlinearities of the system above, we use a fixed-point ap- 
proach. We suppose that, at a given discretized time f" = «Af, 
we know the discretized variables h" = h(x, t"), c" - c(x, f"), 
and u" = u(x, ?")• Let j denote the index referring to the fixed- 
point iteration. Remark that, for j - 0, the generic variable 
equal to a". At the j-th fixed-point iteration, we first 
compute the discretized piezometric head by solving 

Eq.Eli.e., 



Figure 13: Downwards injection at a reference molarity C"' 



0.5 mol/L. 



1 + yc"+^.J V ^' 



= 0, 



(11) 



with the boundary conditions prescribed above. 

Then, using Eq.[8] we compute the discretized advection field 

u"^''-'^' as follows 



u 



n-Hj+l 



K 



1 + yc 



j+i , «-i-i,/„ 
■' + ec ■'e. 



(12) 



Finally, we compute the discretized concentration field c"+''^^' 
using the discretized ADE|3] 



At 



,y«+lj+lj^«+l,i+l^(J3) 



with the boundary conditions above and 
=Q'|u"+'-'+'|/6' + £>o- 



(14) 



Note that we have written the discretized equations includ- 
ing also viscosity variations and molecular diffusion, for sake 
of completeness. In actual calculations, these two terms are 
dropped. The convergence of the fixed-point iterations loop 
is achieved when the error on the generic variable, - 
is below a given threshold. Then, the following time 
iteration is started, f"+' - t" + At. 

The space discretization of the equations is performed via 
the finite element solver available in CAST3M. Within this 



Contaminant profile Cdt) measured at section t = 11 cm, as a function of time 
[h]. Experimental data (averaged over 5 runs) are in blue, model outcomes 
(averaged over 5 initial random flux perturbations) are in red. Error bars for 
both curves correspond to one standard deviation. 



code, several finite element types can be selected to dis- 
cretize the variables and the space-dependent coefficients. 
In this work, both variables and coefficients are approxi- 
mated by_jj^ing_^^agrange complete quadratic polynomials 
(Q2) ( Dhatt and Touzot 1981b . These polynomials involve hn- 



ear combinations of 1, x, y, x^, xy, y , xy , x y, x y terms . 



in the frame of the reference element tohatt and Touzot 198lh . 
For our calculations, we typically used 10 elements along the 
horizontal direction and 160 along the vertical direction. This 
ratio was established on the basis of the aspect ratio of the col- 
umn, i.e., HI2r - 16. 

Finally, the solution of the hnear systems arising from the 
Poisson equation for h (Eq. [TTJ and the discretised ADE equa- 
tion for c (Eq. [T3T l is obtained via the Stabilized Bi-Conjugate 
Gradient method coupled with an ILUTPpreconditioner, whose 
dimension is take n 1.5 times larger than ffie matrix involve d in 
the linear system (ISaad 1996t ISaad and van der Vorst 20001) . 

The time step At is in principle arbitrary, and must be cho- 
sen small enough to attain convergence. In our numerical 
tests, it turns out that this condition is typically ensured when 
At < min (sAx, sAy}, where Ax and Ay are the spacing of the fi- 
nite elements in horizontal and vertical directions, respectively, 
and s ^ 0.02. Tests of convergence were performed and gave 
satisfactory results. 



5. Testing model on experimental data 

In the following, we address the comparison between the ex- 
perimental results described in Sec. |2] and the numerical simu- 
lations of the model, as exposed in Sec.[3]and|4] The measured 
hydraulic conductivity is K ^ 1.59 • 10"^ m/s, with po - 998.3 
Kg/m-' and //q = 0.997 • lO"-' Kg/ms, which allows estimating 
ffie intrinsic permeability k ^ 1.62 • 10 nr. From the data 
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(a) (b) 

Figure 14: Concentration tields corresponding to downwards (left) and upwards 
(right) contaminant injection, at a given time. Stable and unstable fronts appear, 
depending on the density contrast between the injected (displacing) and resident 
(displaced) fluids. Color scale ranges from C/C"'°' = 1 (red) to C/C'""' = 
(blue). 



(a) (b) 

Figure 15: Steam lines corresponding to the cases presented in Fig. 1141 i.e., 
downwards (left) and upwards (right) contaminant injection, at a given time. 
Stream Unes are distorted due to the coupling between advection and concen- 
tration fields. Distortion is evident where the contaminant plume is located 
(cf. Fig.O. 



in Table [T] we have estimated e ^ (Ap/po)/Ac = 0.12. Un- 
certainties on experimental values have been attributed to the 
parameters e and a, which have been adjusted so to improve 
the fit of model outcomes with respect to measured data. 

Since our emphasis is on the spatial and temporal features of 
the contaminant plume dynamics, we commence by consider- 
ing the evolution of concentration profiles measured at several 
heights along the column, as a function of time. In Fig.|8] we 
display as a representative case the results for downwards in- 
jection at C™"' - 0.2 mol/L. As experimental data correspond 
to section-averaged X-rays measures, simulated profiles have 
been as well obtained by averaging the (2d) concentration field 
c(x, f) over the transverse direction x. To this aim, after com- 
puting the field c(x, f) from the model, we define the section- 
averaged normahzed concentration C(j, f) = J c(x,t)dx/B to 
make comparisons with experimental data possible. 

As a general remark, the specific shapes of the unstable fronts 
measured i n experiments ( v yhere the hea vy fluid overlies the 
lighter one ( Wooding 1969t Tavlor 19"5 ^) are not expected to 
be exactly reproduced by transverse-averaged simulated pro- 
files; indeed, the observed profiles intrinsically depend on the 
random small-scale perturbations (such as different grain sizes, 
or preferential flow streams) encountered by the plume at in- 
jection and all along the column. These perturbations, which 
are ubiquitous even in (macroscopically) homogeneous media, 
vary at each experimental run and are ultimately responsible for 
the extension of the unstable front. 

From the point of view of simulations, instabilities 
in the nonlinear system can be triggered by small ran- 
dom perturbations such as the unavoidable numerical 
noise that is produced by integrating the model equa- 



tions ("Tartakovskv and Mea kin 2005 ). In practice, however, 
this procedure could require a time span that is longer than 
the simulation time, since the apparent geometric symmetries 
(the rectangular domain with a flat injecting surface) hinder the 
growth of such perturbations. Therefore, in order to force the 
appearance of interfacial instabilities, it is expedient to perturb 
the injected flux by superposing an additive white noise. In 
our computations, we have chosen a relative noise amplitude of 
1%: these small deviations are sufficient to trigger instabilities 
within the first simulation time steps. Physically, the additive 
noise could be justified on the basis of the small-scale fluctua- 
tions that are generated when injecting the contaminant inside 
the column: the adopted amplitude is well within instrument 
accuracy. 

As discussed in Sec. |2l the stronger the solution molarity 
C™"', the stronger the effects of stochastic fluctuations on con- 
taminant concentration profiles. However, the overall behavior 
at C™°' = 0.2 mol/L, and especially the skewness of the curves, 
is well captured by the model at each position along the column. 
In particular, we may conclude that the approximations that we 
have introduced (such as considering a simplified geometry, ne- 
glecting molecular diffiision and viscosity variations) are a pos- 
teriori justified on the basis of the comparison between model 
and data. In Fig.|9] we perform a similar analysis for upwards 
injection at C'"°' - 0.2 mol/L, for the same step-injection flow 
conditions. Also in this case, a good quantitative agreement is 
found between simulated and measured concentration profiles. 

Up to a normalization factor, the transversally-averaged con- 
centration C(f, t)dt can be interpreted as the distribution func- 
tion of the contaminant plume passage times t{{) at various 
heights y - t along the column. Then, the comparison between 
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model and experimental data is substantiated in Figs.fTOlandfTTI 
by computing the first two cumulants of the passage times, i.e.. 



pollutant plume, even at strong molarities. Similar results have 
been obtained also for upwards injection (not shown here). 



Ki[t({)] = J tC{{,t)dt 

KlW)] = j it-KifC({,t)dt, 



(15) 



for the cases of Figs.[8]and|9l respectively. The quantity ki [t(£)] 
[h] is the mean of the distribution and defines the average time 
required by the contaminant plume to reach height { from the 
injection position; the quantity K2[t{{)] [h^] is the variance and 
defines the (squared) spread around the average time. Model 
outcomes provide accurate estimates of the moments computed 
from experimental data. Note in particular that the quantity 

V - {/Ki[t{i)] identifies the average velocity of the pollutants 
in the porous media: to a first approximation, it turns out that 

V ^ vo, where vo - 6 cm/h is the nominal plume velocity in ab- 
sence of flow-transport couplings. This implies that the average 
velocity is minimally affected by density variations (for the ex- 
perimental conditions considered here). On the contrary, fluctu- 
ations around the mean velocity do not simply average out, and 
are responsible for the increased dispersion in correspondence 
of the unstable front between resident and displacing fluids, i.e., 
the skewed concentration profiles. 

The deviation from Gaussianity of the concentration profiles 
shown in Figs. [8] and |9] is better elucidated by means of the 
third cumulant K^lHi)] - J{t - KifCit, t)dt [h^], which is dis- 
played in Fig. [T2I For Fickian transport, Gaussian (symmetri- 
cal) concentration profiles would lead to vanishing KT\t{()\, on 
the contrary, the moment estimates computed from experimen- 
tal data indicate that downwards injection leads to negatively 
skewed profiles (/C3[f(^)] < 0), whereas upwards injection leads 
to positively skewed profiles (/C3[f(^)] > 0). Remark in particu- 
lar that the absolute value of the third moment is an increasing 
function of column length i, which means that deviations from 
Gaussianity become more apparent as the contaminant plume 
moves through the porous material. Model estimates, which 
are displayed in the same figure, are in good agreement with 
experimental data. 

Finally, in order to further support our analysis, we analyze 
the concentration profiles at the last measure point along the 
column ({ - 77 cm) at C'""' = 0.5 mol/L, for downwards in- 
jection at the same flow conditions as above. We know from 
previous discussions that concentration profiles at this molarity 
show apparent fluctuations that hinder reproducibility; model 
simulations show a similar behavior Then, in order to make 
comparisons possible, we choose to average measured concen- 
tration curves over several experimental runs, and to average 
model outcomes over several random realizations of the initial 
flux perturbations. In Fig. [T3] we display the averaged pro- 
files as well as the error bars (corresponding to one standard 
deviation). Although the model can not capture the precise 
shape of the single concentration profiles, a remarkable quanti- 
tative agreement is found for the averaged realizations. In other 
words, the model is actually capable of correctly predicting the 
average displacement and the spatial extension (spread) of the 



6. Qualitative features of variable-density transport 

The finite elements numerical tool described in Sec. |4] pro- 
vides a practical means of exploring the contaminant dynamics 
inside the column, and can also be used to extract informations 
that are not easily accessible (or even not available at all) from 
experiments. In particular, we are interested in exploring the 
general qualitative features of the interfacial behavior In fact, 
the X-rays measures correspond to averaging the concentration 
profiles in the transverse dimension at any fixed height along 
the column, and thus hide the fine-scale properties of such phe- 
nomena as fingering and instabilities. 

Typical simulated concentration fields c{x,y) correspond- 
ing to variable-density contaminant transport are displayed in 
Fig. [14] at a fixed time t - 2.1 h, for both downwards and 
upwards injection (with duration T - 1.25 h). For this sim- 
ulation we chose C'"°' = 0.3 mol/L, with column geometry 
H - AQ and B - 15. We remark that, because of the flux 
shape (with a finite time duration), two fronts appear at the 
interface between injected and resident fluids: one is unsta- 
ble (owing to the density contrast with respect to the resident 
fluid) and displays fingers, whereas the other is stable. In 
correspondence of the stable front, a standard Gaussian con- 
centration profile is found. The unstable front gives rise to 
an enhanced extension of the mixing (dispersive) region of 
the injected plume, in agreement with experimental observa- 
tions (D'Angelo et al. 2008; Debacq et al. 2001). Note that the 
extension of the fingers is not entirely symmetric upon reversal 
of the flow direction. 

We observe that, once triggered, unstable fingers can grow 
only provided that the contaminant molarity is sufficiently 
strong to overcome the opposing effects of dispersion, which 
acts as a smoothing process on the concentration profiles. When 
transport is Fickian (i.e., at weak molarity), our simulations 
show that dispersion dominates and the small flux perturba- 
tions are rapidly reabsorbed, so that after a few time steps 
the contaminant profiles become perfectly smooth. In par- 
ticular, the shape of the concentration profiles becomes inde- 
pendent of the specific realization of the initial flux perturba- 
tion. Th is is indeed cohere nt with the experimental findings 
of, e.g., dOe Wit et al. 2005h . The number of fingers that can 
actually appear and survive the dispersive smoothing strongly 
depends on the geometry of the column: in general larger diam- 
eters lead to an increased number of such fingers (as naturally 
expected), for a given value of the molarity (Taylor 1950). 

The advection field u(x, t) is also deeply affected by den- 
sity variations. For illustrative purposes, in Fig. [15] we plot the 
stream lines corresponding to the case presented Fig. [141 The 
stream function i/'(x,y) is defined as the scalar function whose 
contour lines are the stream lines. For two-dimensional flows 
satisfying Eq.|2l i^(x,y) is such that 

_ dilj{x,y) 
dy 
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dx 



(16) 



at any fixed time f, Ui and u,. being the horizontal and vertical 
components of the velocity field u(x, y), respectively. Then, it 
follows that Ai/'(x,y) - -V x VL(x,y). By definition, the stream 
lines are instantaneously tangent to the velocity vector. From 
Fig. [15] it is immediately apparent that stream lines for both 
downwards and upwards injection are perturbed where the con- 
taminant plume is located, by virtue of the coupling between 
concentration and advection. This is to be compared with the 
expected behavior of Fickian transport, where the vertical com- 
ponent of the velocity field is constant (and equal to vq) and the 
transversal component vanishes by virtue of symmetries, so that 
the advection field is homogeneous and stream lines are parallel 
and oriented along the vertical direction. 

Finally, observe that the time duration T of injection can 
pla y an important ro l e, as experimental ly investigated, e.g., 
in jwood et al. 2004 lOe Wit et al. 2005b . Indeed, while for 
two semi-infinite fluids only their mutual interface matters, in 
our case two distinct fronts arise. At the unstable front, dense 
fluid fingers move downwards under the influence of grav- 
ity, whereas light fluid fingers move upwards due to buoyancy 
forces (this is particularly evident in Fig. [14] for downwards 
injection). We have numerically verified that, if the time du- 
ration of the injection is short (so that the spatial extension 
of the contaminant plume is limited), light fluid fingers might 
reach the stable front at the opposite end and disrupt the in- 
tegrity of the plume. At longer time scales, dispersion even- 
tually takes_ovei_and_the fingering dies out, because of dilu- 
tion ( jPe Wit et al. 2005h . These phenomena have been tested 
numerically, but do not actually occur for the experimental con- 
ditions under investigation here. 



7. Conclusions 

In this work, we have addressed the problem of variable- 
density transport of a miscible fluid, with the aim of ascer- 
taining the spatial and temporal features of such contaminant 
dynamics in porous media. This topic has been the subject of 
intense research activity, by virtue of its relevance in such tech- 
nological challenges as polluted site remediation, and enhanced 
oil recovery. Our investigation has been supported by novel ex- 
perimental results. 

Data have been collected by means of a vertical column setup 
filled with homogeneous saturated sand, coupled with a dichro- 
matic X-ray source and associate detector. As contaminant 
spills in groundwater typically have a finite extent, our anal- 
ysis has focused on finite-duration pollutants injections. Exper- 
imental results show that the concentration profiles are skewed 
because of gravity effects, and that the sign of the skewness 
depends on the flow direction. 

In order to gain a deeper insight on the underlying physi- 
cal phenomena, we have briefly reviewed the nonlinear equa- 
tions that govern variable-density transport and determined 
their numerical solutions by resorting to the finite elements 
code CAST3M. We have then compared the simulation results 



with the experimental data: the model actually captures the 
salient features of the contaminant dynamics and displays a 
good quantitative agreement with measurements. 

Finally, the finite elements model has been also used as a 
means of exploring the behavior of the interfacial dynamics 
between the dense contaminant plume and the lighter resident 
fluid that saturates the column, whose fine-scale properties are 
not accessible by experiments. Such interfacial dynamics is ul- 
timately responsible for the skewed shape of the contaminant 
profiles within the column. For downwards injection, gravita- 
tional instabilities appear at the front of the plume, and the en- 
hanced mixing is such that contaminant reaches the outlet ear- 
lier than expected when adopting the standard ADE approach. 
At the opposite, for upwards injection the instabilities appear at 
the rear of the plume, and contaminant reaches the outlet of the 
column later than expected. 

In order to perform the comparison between model and data, 
some simplifying assumptions have been introduced, which 
surely deserve further consideration. For instance, it would be 
important to isolate the effects of viscosity (running tests on 
horizontal columns) and molecular diffusi on (at very slow flo w 
rates or with ad hoc devices such as in dKirino et al. 2009h ). 
Moreover, we should separately test the relevance of the 
Boussinesq approximation and extend the numerical simula- 
tions to 3(i cylindrical geometry. 

Many issues remain to be addressed, in view of such applica- 
tions of our laboratory-scale study as to polluted site remedia- 
tion. A fundamental question is whether the present results can 
be extended to larger domains (field-scale measures) via simple 
scaling arguments, or more complex mechanisms and couplings 
must be taken into account. In particular, it has been pointed 
out that t he validi ty of Darcy's and Fick's laws might be ques- 
tionable JHassanizadeh and Leijnse 19951) . As larger scales un- 
avoidably involve complex spatial patterns, preliminary experi- 
mental work will concern the case of heterogeneous and/or non- 
saturated media, where deviations from Fickian behavior due 
to the non-homogeneous structure of the porous material may 
compete with those of density. In this context, preliminary con- 
siderations seem to suggest that a suitable modeling tool could 
be provided by an extension of the conti nuous time random 
walk formalism to nonlinear transport (.Zoia et al. 20091) . 
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